Ground State Properties of Fermi Gases in the Strongly Interacting Regime 
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The ground state energies and pairing gaps in dilute superfluid Fermi gases have now been calcu- 
lated with the quantum Monte Carlo method without detailed knowledge of their wave functions. 
However, such knowledge is essential to predict other properties of these gases such as density ma- 
trices and pair distribution functions. We present a new and simple method to optimize the wave 
functions of quantum fluids using Green's function Monte Carlo method. It is used to calculate 
the pair distribution functions and potential energies of Fermi gases over the entire regime from 
atomic Bardeen-Cooper-Schrieffer superfluid to molecular Bose-Einstein condensation, spanned as 
the interaction strength is varied. PACS: 03.75. Ss, 21.65.+f, 02.70.Ss 
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Recent progress in experimental 
theoretical methods 0, H, 0, E3 has generated great in- 
terest in the properties of dilute Fermi superfluid gases. 
Such gases are also of interest in studies of astrophysical 
objects such as neutron stars [llj and in nuclear physics 
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The Hamiltonian of these gases has the standard form 
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The range of the interatomic potential v(rij) is much 
smaller than the interparticle spacing in the gas, and 
only the s-wave scattering length a of the interaction is 
relevent. Weak attractive interactions have a small neg- 
ative a which increases in magnitude as the interaction 
gets stronger. The a — > — oo as we approach the bound 
molecular state. On further increase of the interaction 
strength, a goes discontinuosly to +oo and then smoothly 
to as the molecule gets more tightly bound. 

Usually, the dimcnsionless quantity l/akp is used to 
characterize the gas. When the interaction is weak and 
attractive, l/akp — > — oo, and we have a BCS super- 
fluid gas with gap A ~ e f/OM (BCS regime). It has 
l/akp << 0, A << the energy per particle Eq/N which 
is positive and less than the Fermi gas energy Efg = 
| w . When the interaction is strong, l/akp >> 0, 
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we have tightly bound molecules with energy E mo i, and 
Eq/N k E mo i/2 « —A. The Bose molecules are con- 
densed in the zero momentum state (BEC regime). In 
the intermediate regime (—0.5 < 1/afcp ^ 0.5) we seem 
to have a smooth transition or crossover from BCS su- 
perfluid to BEC. 

The problem of calculating the ground state energies 
and pairing gaps of superfluid gases has been solved 
with the fixed node Green's function Monte Carlo (FN- 
GFMC) method HUGH- To begin with, we review this 
method and the problem it faces for computing observ- 
ables other than energies. 

In FN-GFMC a trial wave function \I>t(R) is evolved 
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We use R = ri, T2, ■ ■ . ;ri>,r2<, . . ., to denote the con- 
figuration of atoms in the gas, and particles 1,2,... 
have spin up and l',2',... have spin down. The sub- 
script FN denotes that the propagator is constrained 
such that the propagated wave function has the nodal 
surface of Wt(R-) at all r. The energy Et is adjusted 
to keep the norm of the wave function constant. At 
large r, the evolved $(t, R) converges to the lowest en- 
ergy state of the system having the nodes of ^y(R), and 
E T = (*(r)|W|*(r)). Without the fixed node constraint 
it will converge to the exact ground state, but uncon- 
strained fermion Monte Carlo calculations become im- 
practical due to uncontrolled growth in sampling errors. 
This is the well known fermion sign problem. From now 
on, we assume that \1/t(R) is real and that r is large 
enough to approximate the limit r — > oo. 

Following Kalos ^4jj the mixed expectation value 
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is calculated using Monte Carlo sampling techniques. 
Since TL commutes with the evolution operator we have 

(H(T))m " ed = (*(r/2)|*(r/2)) = {U{t/2)) ' (4) 

The (H{T)) m i xe d converges to the energy of the low- 
est enegy state with the nodal surface of \&t(R)- By 
the variational principle it is > the ground state energy 
Eq. The nodal surface of \&r(R) is varied to minimize 
(H(t)) mixed- The minimum gives an accurate estimate 
of Eq provided the variation is general enough. 

This procedure assures that the nodes of ^y(R) are 
near optimum, i.e. close to the nodes of the exact 'I'o(R)- 
However, the >Ft(R) itself can otherwise be very different 
from the \&o(R)- For example, in references || and 
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FIG. 1: The trial, mixed and extrapolated (?Tt( r ) with the 
LOCV ^t(R.) are compared with the trial g-\f(r) with the 
^optimCR-)- Note that the mixed and trial pair distributions 
are the same for ^ op tim(R.) ■ ro is given by ^■nr'^p = 3. 



FIG. 2: Optimized ftt(r) for different values of afcp. 



we use 



*t(R) 



(R) , 
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where $bcs(R) is a generalized Bardeen-Cooper- 
Schrieffer wave function, and its nodes are optimized by 
minimizing {TL{t)) mixed- The /fjXry) is a nodeless pair 
correlation function between spin up and down particles. 
The ^H{t)) mixed does not depend upon the choice of 
ftl( r ij') m the limit t — ► oo. The fn(rij>) is used to 
reach this limit quickly, and to reduce the variance of 
the stochastic evaluation of {TL{T)) m i xe d- Note that the 
commonl y u sed Jastrow pair correlation function is also 
nodeless |l£j, but it acts between all pairs: ||, ft an d 
H, and it is not useful in superfluid gases. 

Mixed expectation values of other observables, 
(C)mixed, can be easily calculated with GFMC, but they 
are more difficult to interprete when [0,Ti] ^ 0. If one 
assumes that |^(t)) = \^t) + |<^), then the desired ex- 
pectation value 
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Here, the trial estimate (O) trial = (^t\0\^t) / {^t\^t}- 
When 5^ is small, the extrapolation 



(0(t)) extrap. ~ 2{0(t)) mixed ~ (O)trial 



(7) 



can be used to estimate (O). 

However, in the strongly interacting regime the 5^ is 
not necessarily small, and this extrapolation may not 
be valid. For example, we consider the pair distribu- 
tion function (r) between parallel spin particles in the 
akF — > ±co limit. The mixed, trial and extrapolated 
values of gn(r) obtained from the 'I't(R-) of Ref. Q and 
are shown in Fig. ^ At small r, the extrapolated 
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FIG. 3: Optimized f n (r) (dashed line) and LOCV / u (r) 
(continuous line) for different values of akF. 



.9Tt( r ) < indicating invalidity. These and all the other 
results presented in this work are obtained from Monte 
Carlo computations using 14 particles in a cubic periodic 
box. As discussed in Ref. Q and , a periodic box with 
14 particles provides a fair approximation to the uniform 
gas. The cosh u(ry/) with (jltq = 12 is used 9] to approx- 
imate the interaction between spin j j pairs. 

In principle, the pair correlation functions, /|j.(r) and 
fll(r) = fn(r) in ^j'(R) can be obtained by minimiz- 



TABLE I: Summary of the results in units of Efg 
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{"^) optim 


-1 


0.792(4) 


0.818(3) 


0.808(4) 


-0.55(3) 


-0.54(3) 


-3 


0.635(6) 


0.85(3) 


0.70(2) 


-1.8(1) 


-2.0(1) 


-10 


0.494(7) 


0.68(3) 


0.53(1) 


-3.5(2) 


-3.2(1) 


oo 


0.414(5) 


0.62(3) 


0.46(1) 


-3.9(1) 


-4.0(2) 


10 


0.32(1) 


0.57(6) 


0.39(1) 


-4.8(1) 


-5.0(1) 
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-0.00(1) 


0.4(1) 


0.11(3) 


-7.0(1) 


-7.3(3) 
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-0.34(2) 


0.2(1) 


-0.18(3) 
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-9.2(3) 
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-2.37(3) 


-0.1(1) 


-2.01(3) 


-19.0(4) 


-18.0(6) 



3 



ing the trial energy ft) trial- However, this variational 
problem has been approximately treated in most quan- 
tum Monte Carlo calculations. In Ref. 0,0 a simple and 
crude method called LOCV, based on constrained mini- 
mization of the leading two-body cluster contribution to 
(ft) trial [H is used - In this method, f^{r) = fa{r) = 1, 
and fn(r) satisfies the two-body Schrodinger equation 

- ^V 2 / u (0 +v(r)f U (r) = Xf n (r) , (8) 

at r < d. The boundary conditions are: /fj_(r > d) = 1 
and fi\(r = d) = 0. The healing distance d serves as the 
variational parameter. The trial energies obtained with 
variational Monte Carlo (VMC) calculations using the 
optimum healing distance d are compared with the FN- 
GFMC (ft) mixed in Table Both calculations use the 
optimum ^sc^R) found by minimizing the ft) mixed 
in Ref. 0- The trial energies are well above ft) mixed, 
particularly at 1/akp > 0. This shows that the LOCV 
pair correlation functions are far from those in the exact 
i n the strongly interacting regime. 
Here we present a new and simple method to optimize 
the pair correlation functions in the trial wave function 
using GFMC results. The optimized trial wave functions, 
denoted by \& op ii m (R) are presumably close enough to 
^o(R-) so that is small and Eq. [7| provides a fair ap- 
proximation. The method can be improved and 5^ can 
be further decreased by including higher order correla- 
tions corresponding to triplet, quadruplet, etc. However, 
in the present work we consider only pair correlations for 

^ optim (R) ■ 

The trial pair distribution functions can be expressed 
as [13 

{g x {r))trial = fx(r) t x (rJ n ,f n ,$ BC s(R),p) , (9) 

where x can be tt or Tl and tx is a complicated function 
of r, /tTj/uj ^BCs(R) and gas density p. It is diffi- 
cult to calculate it exactly except by numerical methods. 
However, t x (r) contains many-body integrals, and is a 
relatively smooth function of r. 

Our method to optimize f x (r) using GFMC is iterative. 
Let (gi n \r)) m ixed and {g x n \r)) trial be obtained from the 
n-th trial fj^ (r) using the optimum &bcs which docs 
not depend on f x (r). We start with the LOCV ap- 
proximation providing the fx(r), but one could start 
with any other choice of fx and converge to the same 
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(R) . The next improved fx (r) is chosen as 
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If the difference between /J^(r) and f x 2 \r) is small, we 
can assume that the t x (r) functions do not change much. 



In this case {gi 2 \r)) tr iai ~ {9x\r)) mixed- Otherwise, by 
iterating this process one easily converges to an f x n \r) 



such that 



(11) 



Usually, the convergence within statistical errors can be 
reached within 3^4 iterations and it doesn't seem to 
depend on the strength of the interaction. In practice, 
{g x (r)) mixed and (g x (r)) trial have Monte Carlo sampling 
errors. We approximate the square root of their ratio 
(Eq. I10|) by a smooth function of r chosen as cos(pir + 
p 2 )e~ r / p3 + l, and vary the parameters P1-3 to best fit the 
Monte Carlo values. One iteration step typically takes 
about 10 hours in a Pentium 3.0 GHz based workstation. 

The VMC energies with the ptim(R-) are much closer 
to the FN-GFMC energies (Table G}. In principle, the 
optimization of f x {r) should have no effect on the FN- 
GFMC ft)mixed; in practice the ft) m ixed seems to get 
lowered by ~ 2 ± 1 % after optimization presumably 
because the limit r — » oo is easier to reach with the 
^ optim(R) ■ The effects of the optimization are also 
seen in the reduced error bars of the energy estimates: 
5 ft) 



optim 



< 5 ft) trial (Table Pi for the same number of 
Monte Carlo samples. In addition, the Et (Eq. [5J, which 
typically has larger fluctuations, becomes indistinguish- 
able from ft) mixed- 

The pair distribution functions (g x {r)) mixed are de- 
termined by the many body probability distribution 
given by ^ optim (R) v I / (T, R), while the (g x {r))o P tim are 
for |* optlm (R)| 2 . Note that * opiim (R)*(r, R) « 



^ optim (R)*o(R) > since the nodes of \I/ p t i m (R) have 
been varied to match those of \&o(R)- 

Extending the above method, if we can match the 
mixed and optimized trial distributions for all, pair, 
triplet, quadruplet, . . . distribution functions, then we 
can assume that ^ optim (R) — 'I'o(R). However, here 
we approximate the exact \&o(R) by \E , op ti m (R) using 
^scs(R) and pair correlation functions only 



TJ fH U m (r i/f )$ BCS (R) - (12) 
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The validity of Eq. II II ensures that the present optimiza- 



tion method will converge to / : 



O) 



f° pUm (r) and thus 



^ojrfim(R) is as close to \I/o(R) as its form (Eq. I12f> al- 
lows. If higher order correlations have negligible effects 
on the wave functuion, we should expect 
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optim 
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E 



(13) 



However, in the interesting regime of akp ~ oo, Table 
[J] shows that the ft) optim is larger than the ft) mixed 
by ~ 10 %. This suggests that the form of the present 



4 



1000 




°\) 0.5 1 1.5 

r/r 
' o 

FIG. 4: Comparison of the pair distribution function <?tl( r ) 
with the radial probability distribution of the molecule in the 
a > regime. 

^'optim(R-) is not sufficiently general. An improved ap- 
proximation could be obtained by including products of 
triplet correlations Fp(rij ,rjk,rki) for ||| and J. J, J., and 
FAi(rij,r jk > ,rk'i) for TTI and ||| triplets in the wave 
function. We believe that the present method can be 
generalize to determine the optimal forms of three-body 
correlations by making 

93,x(Tij i Tjk !^ki)mixed — QS.xij'ij : Tjk > Tki)optim : (1^0 

where g% denotes three particle distribution functions. 
The true ^ can also have backflow correlations [T^ : 
however, they change the nodal surface and have to be 
optimized by minimizing (H) mixed- 

The main difference between the optimum pair corre- 
lations and those of Ref. @ and is in f-\i(r). In LOCV 
we have /n( r ) = E because in two-body clusters there 
is no interaction between parallel spin particles in dilute 
Fermi gases. However, many body effects generate an ef- 
fective repulsion between parallel spin particles and the 
optimum ff[(r) is < 1 at r < 1.5r as shown in Fig. |2] 

The optimum and LOCV are generated by the 

strong two body attraction in f J, pairs and have quali- 
tatively similar shapes (Fig. |3J). For 1/akp « —1, the 
LOCV fn(r) is near optimum. For stronger interactions, 
it is larger than the optimum at r ~ (Fig. |3J. 

The expectation value of the potential energy, (V) = 
(J2i j> v ( r ij')) can easily be calculated from the GFMC 
(mixed) and VMC distributions using 'f optim- The calcu- 
lated values of the potential energy are given in the last 
two columns of Table [I] Apart from statistical fluctua- 
tions, the mixed and the optimum pair distribution func- 
tions are the same. Therefore no extrapolation, such as in 
Eq. [7| is necessary for calculating (V) using ty ptim(R-)- 

Only when 1/akp > 0, we can have bound states 
with normalized radial wave functions R(r). We de- 



fine R'{r) — yJ~R(r) so that <jw( r u) = R' 2 { r ) is nor- 
malized analogous to g-\i{r), and the two are compared 
in Fig. 0] (six(r) is normalized such that g-\i(r) — > 1 
for r — > 00). When 1/akp — > + , we know that 
9moi( r ) = R' 2 (r) — » (infinite pair size), but gn(r) > 1. 
So g m oi( r ) and g-\i(r) are qualitatively different when a 
is large. However, when the interaction is stronger and a 
becomes positive and small we expect a gas of molecules 
in which <7fj.(r) ~ gmoi( r ) at r < the size of the molecule. 
Fig. 01 shows that the superfluid may well be approxi- 
mated by a gas of molecules with BEC for 1/akp > 1/3. 
The molecule size for 1/akp = 1/3, is 1.21ro and for 
1/akp = 1, it is 0.38ro. Due to the many body effects 
in g^i(r), it is meaningless to compare beyond these dis- 
tances. In fact, for all akp > 0, the g m oi( r ) — * while 

9u( r ) - > 1 as r - > °°- 

In conclusion, the proposed method allows us to op- 
timize separately the BCS and the pair correlations in 
dilute Fermi gases. The BCS and fn(r) correlations 
are most important, however in the strong interaction 
regime, the /n( r ) can not be neglected. 

The studies of momentum distributions and density 
matrices of the superfluid gas may now be possible using 
the optimum ^r(R), and are in progress. 
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